#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <pthread.h>
#include "rand_bm.h"

// N_THREADS e' o numero total de threads que serao executadas
// REAL_THREADS e' o numero de threads menos um, porque eu vou aproveitar o main como uma thread
#define N_THREADS 4
#define REAL_THREADS (N_THREADS-1)

double *trials;

int M;		// M : número de iterações
double S;	// S: valor da ação
double E;	// E: preço de exercício da opção
double r;	// r: taxa de juros livre de risco (SELIC)
double vol;	// σ: volatilidade da ação
double T;	// T : tempo de validade da opção
int range;

int mean_sum[N_THREADS];	// somas parciais

// retorna o maior entre a ou b
double max(double a, double b) {
	if(a > b) return a;
	return b;
}

// calcula a media aritimetica de trials, mas cada thread ja fez uma soma parcial 
// "mean_sum", entao e' soh somar as parciais e dividar pelo numero de elementos
double fmean() {
	double sum = 0;
	int i = N_THREADS;
	while(i--)
		sum += mean_sum[i];
	return sum/M;
}

// calcula o desvio padrao de trials, dado uma media mean
double fstddev(double mean) {
	double aux;
	double sum = 0;

	int i = M;
	while(i--) {
		aux = trials[i] - mean;
		sum += aux * aux;
	}
	return sqrt(sum/M);
}

// calculo do blackscholes
void* blackscholes(void *valor) {	// valor e' a posicao inicial que a thread vai trabalhar no vetor trials.
	struct BoxMullerState state;
	initBoxMullerState(&state);

	double t;
	int b =((int)valor);	// b assume o valor da posicao inicial, so' fiz isso para evitar parenteses no codigo
	int index = b/range;	// index serve para colocar as somas parciais na posicao certa do vetor mean_sum
	mean_sum[index] = 0;	// somas parciais do trials calculados

	int i = b;
	for(; i<b+range; i++) {
//		t = S * exp( (r - 0.5*(vol*vol)) * T + vol * sqrt(T) * (rand() / (double)((1LL<<31)-1)) );
		t = S * exp( (r - 0.5*(vol*vol)) * T + vol * sqrt(T) * boxMullerRandom(&state) );
		trials[i] = exp(-r * T) * max(t-E, 0);
		mean_sum[index] += trials[i];		// ja calcula somas parciais do trial, para usar depois no calculo da media
	}
}

int main(int argc, char *argv[]) {
	scanf("%lf%lf%lf%lf%lf%d", &S, &E, &r, &vol, &T, &M);	// inicializa as variaveis a partir dos dados entrados

	trials = (double*)malloc(sizeof(double)*M);

	srand(time(NULL));

	range = M/N_THREADS;		// range e' quantos elementos de trials, cada thread vai ficar responsavel
	pthread_t t[REAL_THREADS];	// cria um total de REAL_THREADS threads

	int i = REAL_THREADS;		// eu procuro inicializar variaveis e utiliza-las logo em seguida para aproveitar o principio da localidade
	while(i--)			// loops deste jeito sao mais rapidos do que for(i=0; i<valor; i++) porque comparacao com zero nao precisa colocar o valor comparado no registrador
		pthread_create(&t[i], NULL, blackscholes, (void*)(range*i) );
	blackscholes((void*)(range*REAL_THREADS));	// o main executa um blackscholes assim como as threads, em vez de ficar esperando

	for(i=REAL_THREADS; i--; )		// espera as threads terminarem
		pthread_join(t[i], NULL);

	double mean = fmean();		// eu retirei os parametros das funcoes, para nao perder tempo colocando as variaveis na stack, ja que elas sao globais
	double stddev = fstddev(mean);
	double confwidth = 1.96 * stddev / sqrt(M);
	double confmin = mean - confwidth;
	double confmax = mean + confwidth;

	free(trials);

	// imprime o resultado na saida padrao
	printf("S\t%d\nE\t%d\nr\t%d\nvol\t%d\nT\t%d\nM\t%d\n", (int) S, (int) E, (int) r, (int) vol, (int) T, (int) M);
	printf("Confidence Interval: (%f, %f)\n", confmin, confmax);

	return 0;
}

